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Neutrino oscillations in the early universe differ from e.g. solar oscillations in two 
important aspects. First, one cannot neglect neutrino annihilation or scattering in 
the medium for a rather large and physically interesting range of parameters (mixing 
angle, sin 26', and mass squared difference, Sm^). These processes break the coherence 
of neutrino propagation and that is why considerations in terms of wave functions 
become impossible and one has to use the density matrix formalism [|l], |^. Kinetic 
equations for the density matrix of oscillating neutrinos with the account of the so 
called second order effects (proportional to the second power of the Fermi coupling 
constant, Gp) were derived in refs. 0, 0) H- Second, neutrino oscillations in the 
primeval plasma could modify the plasma properties, in particular the refraction 
index, and this in turn would influence the oscillations, so that the problem becomes 
highly nonlinear. The refraction index of oscillating neutrinos in the cosmic plasma 
was calculated in ref. [Q. An important feature of the refraction index is that it 
contains terms proportional to the charge asymmetry of the cosmic plasma. The 
effective potential of a standard neutrino of flavor a can be written as 

.GlT^E ,^, 



V^^ff = ±CitiGfT' + C2"- 



a 
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where E is the neutrino energy, T is the temperature of the plasma, Gp = 1.166 ■ 10 
GeV~^ is the Fermi coupling constant, a = 1/137 is the fine structure constant, and 
the signs "±" refer to anti-neutrinos and neutrinos respectively (this choice of sign 
describes the helicity state, negative for u and positive for z/). According to ref. |Q 
the coefficients Gj are: Ci ~ 0.95, C| ~ 0.61 and G2'^ ~ 0.17. These values are true 
in the case of thermal equilibrium, and otherwise these coefficients are some integrals 
over the distribution functions. In ref. p the coefficient Ci was calculated using the 
present value of the asymmetry 77, which differs from its value in the early Universe 
(before e~^e~ annihilation increased the photon temperature) by a factor 11/4. In our 
calculations we took Gi = 0.345. The charge asymmetry, 77, is defined as the ratio of 
the difference between particle-antiparticle number densities to the number density 
of photons. The individual contributions to 77 from different particle species are the 



following: 

7] = 2r]^^ + 7]^^ + rjy^ +r]e- r]n/2 (for z/g) , (2) 

V = "^Vu^ + Vu, + Vur - ^n/2 (for z/^) , (3) 

and rj for z/^ is obtained from eq. (|^) by the interchange ^ ^^ t. 

The different magnitude of the refraction indices for neutrinos and anti-neutrinos 
may result in more favorable conditions for Va ~^ ^s oscillations compared with z/^ -^ 
Us oscillation, where z/^ is an active neutrino, a = e, /i, r, and Ug is a sterile one. Since 
more z/^ than z/q would be transformed into sterile ones, the lepton asymmetry in 
the sector of active neutrinos would rise and this would further amplify the process. 
The possibility of such instability was noticed in ref. 0, but there it was found, 
on the basis of simplified considerations, that the rise is stabilized when non-linear 
terms in the refraction index become non-negligible, and thus it was concluded, that 
"no large chemical potential will be generated in any point of the parameter space" . 
A similar statement of a small asymmetry was made in ref. M. This conclusion 



was reconsidered in ref. P| (see also refs. M, |TD|) where it was argued that a very 



large asymmetry, even close to 1, may be generated by the oscillations. An even 



more striking statement was made in ref. |TT|, that the asymmetry may not only 
be large, but may have a chaotic sign, so that even domains with different signs 
of the asymmetry may be formed |1^. Similar results were obtained in the recent 



paper |Q, namely, that the asymmetry may reach large values and in some ranges 
of parameters its sign may be chaotic, while in ref. [jl4| chaoticity was not observed. 



However, these results were derived in a simplified way after averaging some essential 
quantities over neutrino momenta or through a solution of momentum dependent 
but simplified equations. At the same time, in "brute force" numerical calculations 
applied to this problem, it is very difficult to distinguish between the real effect and 
a computational instability. Technically there is an essential difference between exact 
momentum dependent calculations and the momentum averaged ones. In the last case 
one has to solve a set of ordinary differential equations, while in the former case the 



corresponding equations are integro-differential ones (see e.g. ref. [^, where exact 
kinetic equations were accurately solved for non-oscillating neutrinos). Momentum 



dependent exact numerical calculations were done in a series of papers ITBI . It was 
shown there that the asymmetry may rise by several (3-4) orders of magnitude but 
still remains small in the range of parameters \5m?\ < 10^^ eV^, assuming a small 
initial asymmetry, 10^^°. Outside this range of parameters the calculations became 
unstable and no definite result was obtained. We extended the domain of stable 
computation to a somewhat larger range of parameters ]T7]: 

\5m^\< 10-^ eV^ , (4) 

and have also not found any large generation of charge asymmetry. In the coinciding 
range of parameters our results are in a reasonable agreement with those of ref. [|16] . 



The attempts to extend the range (^) maintaining the stability of the computational 
process demanded a huge increase of computer time, because otherwise the results 
were chaotic in sign and showed a quickly rising asymmetry. Thus it seems that 
further attempts to extend stable numerical calculations to a wider parameter range, 
in which the huge amplification of the asymmetry was found, will be fruitless and one 
should try to transform the kinetic equations for the oscillating neutrinos analytically 
to such a form that will permit a numerical solution. In what follows we have achieved 
this goal and reduced the problem to the solution of an (almost) ordinary differential 
equation for the neutrino charge asymmetry, which is easy to solve numerically. 

Before presenting the actual calculations we will briefly describe our procedure 
so that it will be easier to follow it. We start from the usual equations for the 
evolution of the neutrino (as well as anti-neutrino) density matrix in a cosmological 
environment. We twice introduce new variables, first x and y given by eq. (|ll|) so 
that the evolution operator in the l.h.s. of the kinetic equations would depend only 
on one variable x, and second, we rewrite all the equations in terms of r given by 
eq. (p2D. It is a natural variable for the description of the behavior of the density 



matrix near the MSW resonance, where it takes the value r = 1, in the limit of 



a vanishing contribution of charge asymmetry to the refraction index. As we see 
in what follows almost all coefficient functions in the kinetic equations depend only 
on r except for the charge asymmetric potential V in eq. (E^) that depends on two 



variables r and y. Using the variable r permits to factor out the large parameter Q 
in eq. (^) related to a large frequency of oscillations. This is true for neutrino mass 
differences above 10~^ eV^. Quick oscillations make numerical computations very 
difficult. Fortunately one may analytically separate the quickly varying functions and 
make an expansion in terms of 1/Q. An essential technical step is to consider charge 
symmetric and antisymmetric elements of the density matrix, p ± p. Antisymmetric 
combination directly enters the refraction index (see eqs. (p],pOD) and working with 
symmetric and antisymmetric functions permits to derive a closed equation for the 
evolution of the asymmetry. As a first step we formally solve equations ([l^ ) -([T7| ) for 
the antisymmetric elements of the density matrix in terms of unknown symmetric 
functions and the integrated charge asymmetry Z in eq. (|20|) . In the limit of large Q 
the corresponding differential equations allow a simple algebraic solution. As a second 
step we substitute the obtained expressions into the charge symmetric equations. For 
the latter we find eigenfunctions that are formal solutions of these equations in the 
case of constant coefficients. However since the latter are not constant we obtain a 
system of differential equations for the coefficient functions in the expansion of the 
solution in terms of the eigenvectors. The equations for these coefficient functions 
are quite simple and can be solved numerically and analytically (both approaches 
give very close results). After that we substitute the found solutions, which contain 
an unknown charge asymmetry Z, back into the antisymmetric equations and after 
integration over momentum of both sides of the equation for (p'^^ — p'^^) we obtain a 
closed differential equation for the charge asymmetry Z. The latter is a function of a 
single variable q = yr and it can be relatively easily solved numerically. Now we will 
describe the same in more detail. 



The basic equations governing evolution of density matrix are: 

i{dt- Hpdp)paa = Fo{psa- Pas)/'2-iTo{paa- feq) , (5) 

i{dt - Hpdp)pss = -Fo{psa - pas)/2 , (6) 

i{dt- Hpdp)pas = WoPas + Fo{pss- Paa)/2-iTiPas , (7) 

i{dt - Hpdp)psa = -Wopsa - Fo{pss - Paa)/2 - iTipsa , (8) 

where a and s mean "active" and "sterile" respectively, Fq = 5m? sin 29 /2E, Wq = 



5m COS26 /2E + K^j, H = J^nptot/'iM'^ is the Hubble parameter, p is the neutrino 
momentum, and feq is the equilibrium Fermi distribution function: 



/e,= [exp(E/T) + l]-^ . (9) 

More precisely instead of the equilibrium function f^q one should use the one with 
a non-zero chemical potential, because scattering and annihilation processes do not 
change lepton number. However if the asymmetry is not large the difference between 
them is not important for our calculations. 

The anti-neutrino density matrix satisfies the similar set of equations with the 
opposite sign of the antisymmetric term in V^tt and with a slight difference in damping 
factors that is proportional to the lepton asymmetry. 

Equations (^^ account exactly for the first order terms described by the refraction 
index, while the second order terms describing coherence breaking are approximately 
modeled by the damping coefficients Vj. The latter are equal to |T8|: 

To = 2ri = Qa ^^^ GpT p . (10) 

In general the coefficient ga{p) is a momentum-dependent function, but in the approx- 
imation of neglecting [1 — /] factors in the electro-weak collision rates it becomes a 



constant [T^ that corresponds to g^^ ~ 4 and Qu ^^^ ~ 2.9 |20|. In the following we will 
use a more accurate value of Qa-, which comes from the thermal average of the complete 
electro- weak rates (with factors [1 — /] included), which we calculated numerically 
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from our Standard Model code |]T5|. This gives us g^^ — 3.56 and g^ ,^^ — 2.5. The 
indices sub-0 are here prescribed to the coefficient functions to distinguish them from 
the similar ones after dividing by Hx (see below). 
It is convenient to introduce new variables: 

X = moRit) and y = pR(t) , (11) 

where R{t) is the cosmological scale factor so that H = R/R and mo is an arbitrary 
mass (just normalization), we choose mo = 1 MeV. In the approximation that we will 
work, we assume that T = —HT, so that we can take R = 1/T. In terms of these 
variables the differential operator {dt — Hpdp) transforms to Hxdx- We will normalize 
the density matrix element^ to the equilibrium distribution: 

Paa = feq{y)[l + a{x,y)], Pss = feq{y)[^ + s{x,y)] , (12) 

pas = pla = feq{y)[h{x,y) + il{x,y)] , (13) 

and the neutrino mass difference Sm'^ to eV^. 

As the next step we will take the sum and difference of eqs. (@)-(|) for u and u. 
The corresponding equations have the following form: 

4 = Fl±, (14) 

a± = -Fl± - 27+a± - 2-f^a^ , (15) 

h'^ = Ul± - VZl^ - 7+/i± - 7_/i^ , (16) 

l'± = ^(a±-s±)-[//i± + VZ/i^-7+/±-7_/^ , (17) 

where a± = {a± a)/ 2 etc, and the prime means differentiation with respect to x. We 

have used W = U ± VZ, 7 = Ti/Hx, and 7-1- = (7 ± 7)/2, where 7_ parameterizes 

^ Other authors find it convenient to express this density matrix formahsni in terms of Paufi 
matrices and a polarization vector, P = (P^, Py, Pz), such that: 

Po 



P- 2 



\ + P -a 



in such a way that PqPz = 1 means that aU the neutrinos are v^, and we have PqPx = feqh, 

PoPy = -feql, PqPz = feq{a ~ s) and Po = feq{2 + tt + s) . 



the difference of interaction rates between neutrino and anti-neutrinos, which is pro- 
portional to the neutrino asymmetry. With the approximation]] ptoj ~ 10.757r^T^/30, 
the expressions for U, V, and Z become: 



2 

U = 1.12 • 10^ COS 2^5m2— + 26.2^, (18) 

y x^ 

V - ^. (19) 

X^ 



Z = W'^'U-J -^y'fe,a\ , (20) 

where rjo is the asymmetry of the other particle species (see eqs. (@,@)) normalized 
in the same way as the neutrino asymmetry (the second term in (pO|)). Here we have 
imphcitly assumed that z/q = i/g- 

We can use total leptonic charge conservation to determine 7_, but as we see 
in what follows, the 7_-terms are either sub-dominant or not important, so we do 
not need a concrete expression for 7_. The most unpleasant contribution, which 
makes it so difficult to solve the symmetric equations numerically, comes from the 
term containing an integral over momentum of the difference {paa — paa) with a large 
coefficient. If the lepton asymmetry is sufficiently small, rj < 10^^, this term can 
be neglected in the symmetric equations, but when it is large its back-reaction on 
the rise of the asymmetry is quite important. All other (asymmetric)^ terms in the 
symmetric equations, e.g. 7_a_, can always be neglected. The only essential terms 
are proportional to V Z because they enter with a numerically large coefficient (see 
below) . 

Let us introduce some more notations. Since the asymmetric term a_ enters the 
expression (|20| ) with the coefficient 10^° we introduce capital letters for the renormal- 
ized asymmetric functions: 

S = 10^°s_, A = 10^°a_, H = 10^°/i_, and L = 10^°L . (21) 



^ This approximation is valid for high temperatures T > 1 MeV. In our case for Sm^ = — 1 eV^ 
and any smaU sin 20 <C 1 we deal with temperatures above 10 MeV for all essential momenta. And 
only for |(5m^| < 10~^ eV^ one will need to take into account a more accurate expression for ptot- 



We will also introduce a new variable: 

r = exVy , (22) 



where ^ ~ 6.5 • 10'^ J\5rn?\ cos 20 so that U vanishes at r = 1 or, in other words, 
the MSW resonance takes place at r = 1 if 5m? < and if the contribution from 
the asymmetric part, VZ, can be neglected. We will divide everything by the factor 
M = 1.12 • 10^ cos 20 ISm'^l x'^/y, so that the coefficient functions now become: 

F = -tan20f« -sin20, U=1/t^-1, j = 6/t'^ , (23) 

where 6 ~ 1/135 is a small coefficient. In what follows we will often use the notation 
7 = 7+. The asymmetric potential, V/M, in terms of these variables has the form: 

\/ = 3.3 ■ 10"^(cos 26 5m^)-^''^yq-^/^ . (24) 

where we have introduced q = yr. The equations for the asymmetric functions can 
now be written as: 

S'/Q = FL , (25) 

A'/Q = -FL - 27A - 2 ■ 10^%_a+ , (26) 

H'/Q = UL--fH- 10^^VZl+ - 10^V^+ , (27) 

L'/Q = ^{A-S)-UH--fL + 10^WZh+-10^^-f.l+ , (28) 

where prime now means derivative with respect to r and: 



Q^5.6-lO^^|(5m2cos20| . (29) 

Due to conservation of leptonic charge the integrated contribution of the last two 
terms in the r.h.s. of eq. (|26| ) vanishes: 



dyy^feM [lA + 10iVa+) = . (30) 

In the equations for H' and LI we neglect the terms proportional to 10^°7_ ~ Z as 
well as Fi^A — S)/2 because they are small in comparison with lO^^VZ ~ lO'^Z. 
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We will solve these equations in the limit of large Q ^ 1, the corrections generally 
being of the order of 1/Q. A formal solution of equations (^) and p8|) (in terms of 
unknown functions Z, /i+, and /+) is: 



L 
H 



^2 + f/2 






(31) 
(32) 



72 + ^/2 

This approximation works ii Z{q) does not decrease too fast with increasing q and is 
even better justified if Z{q) is a rising function of q. 

These solutions can now be inserted into the set of equations for the symmetric 
function, which can be written in the matrix form as V = M.V\ 



V i'+ J 



Q 



F \ 



000 

-27 -F 

-7 f7 

\ -F/2 F/2 -ii -7 y 



( s^\ 



V 1+ J 



(33) 



Here we have used: 



U = U{1- D^a^) and 7 = 7 (1 + D^a^) , 



(34) 



with L)2 = ^2y2 g^^^ ^2 ^ ^2 ^ ^2_ 

We will solve this system of equations expanding the solution in terms of the 
eigenvectors of the corresponding matrix Ai in the r.h.s. of eq. (|33D with the coeffi- 
cients bj that are not constants (as they would be if Ai was a constant matrix), but 
functions of r. For the functions bj we will obtain a set of differential equation that 
can easily be solved numerically and even analytically in the limit of a small sin 26. 
The matrix Ai has 4 eigenvectors with the eigenvalues: 

F27 



/ii 

fJ'2 



2^2 

-27 + 



F2(27-7) 



2[(27-7)2 + f/2] 
-7±if/ 



(35) 

(36) 
(37) 
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where 5"^ = 7^ + t/^. The correction to /i2 is of the order of F"^ when D^ — a^ 7^ 0. 
When the latter quantity is close to zero, the correction may be of the order of F. 
The matrix elements of the symmetric density matrix can be expressed through 



the four new functions hj{T) as: 



s+ 



-&o + hi 
-F 



4a2 



a^a^ 



7 (62 cos fi — 63 sin Q) — If (62 sin i7 + 63 cos Q) 



a+ 



h. 



-bo r^^z:^ ^oi- 



Aa'^a 



?r2 



a' 



-— [7 (62 cos fi — &3 sin Vl) + U (62 sin fi + 63 cos Vt)] , 



cr" 



I, 



bo 
bo 



FUja^ - D^ 
F7(a2 + D^) 



bi TT^ H o (,"2 sm i 2 + 03 cos i 2) 



2a2 



a^ 



F^ a^ - ^2 
^iTT^ "I 9 — (62 COS i7 — 63 sin fi) , 



a" 



(38) 

(39) 
(40) 
(41) 



2a2a2 "2^2 

where Vt' = QU . To the leading order in the small parameter F the function bo 
satisfies the equation^ 



b' 



2^2 



-6n. 



(42) 



We usually neglect terms of the order of F"^, but the one above contains the large 
factor Q and is therefore taken into account. 

The function 61 is small, 61 ~ F"^, and can be neglected. The functions (62,3"^) 
satisfy the equations: 

bo 



(63^)' = -Q7(63^) 



bo 



{ail))' smVl + {(3(j)\ cosil 
(a^)'cosr2— {(3(j)\ sinfi 



(43) 
(44) 



where V' = (a^ - D'^)/a'^, (j) = {a'^ + D^)/(x'^, & = FU/a"^, and /3 = F-f/a^. The initial 

conditions are 6o(0) = —1 and &i, 2,3(0) = 0. 

^Note that this equation can be written as the evolution equation for the sterile neutrino (and 
anti- neutrino) density in momentum space as shown in ref. |9[| (eq. (93)) or ref. ||l9[| (eq. (72)). 
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When D = it is straightforward to solve for the 6-functions numerically. From 
fig. ID it is clear that &i,2,3 are very small except near the resonance. For momentum 
y = 1 and parameters sm26 = 10~^ and 6771"^ = —1, we see that the function 69 
follows the curve exp [— k ■ (arctan(2(g — 1)/S) + 7r/2)] to a high accuracy (remember 
that q = yr). Here k, goes like sin^2^ for small mixing angles. 

The last two equations ( ^^if) can be solved as: 



h^p = -]- r dne~^'^^^^^'^^'\{Ti) (aV^);cosfii- (;30)\infii 
where Fg = Q7 and sub-1 means that the argument of the corresponding function 



(45) 
(46) 



is Ti. When substituting these results into eqs. (^,0) and integrating by parts we 
obtain: 






Q 

2 Jo 



drie 60 (n) 



K 



F sin AQ - — ^ f/5i0i sin AQ - oiipi cos AQ 



Qbo 



2 Jo 



b' 



F cos AQ - — ^ f ai^i sin AQ + /5i0i cos AQ 



Qbo 



(47) 

m 



where AF3 = T^{t) — F3(ri) and AVt = VL[ti) — VL[t) . The integrals can be taken in 
the limit of large Q according to: 

1 $(r) 



dTi^[Ti)e 



-AT-iAn 



Q^~iU 



(49) 



This result permits us to express the function L from eq. (|3lD algebraically through 
the lepton asymmetry Z: 



-10 



10 



FVZjU 



F\a' + D' 
40-25-2 



(50) 



Substituting this result into eq. (p6|) and integrating it with dyy'^ feq{y) / 4:7i'^ we finally 
derive the following equation governing the evolution of the lepton asymmetry Z{q): 



1 dZ 
Z dq 



-5Bq^/^ / dt 
Jo 



t\t' - l)UtqMl/t) / _ ^^a' + D' 



fj^cr^ 



cr^O"^ 



(51) 
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where: 

B = 4.71.10''(cos20|*m=|)'''»[^pf] , (52) 

B, . 2.,.10-^[^)\ (53) 

and where we have introduced the new integration variable t = 1/r, which is pro- 
portional to the neutrino momentum. If we neglect the term proportional to Bi our 
evolution equation (|5T|) coincides with the main contribution of the "static approx- 
imation" in refs. P, ^ (see for instance eq. (65) in [|l9l). However the last term in 



eq. (|5TD stops the rise of the asymmetry Z earlier and at much lower values. 

In the limit of a small Z we may neglect D in the r.h.s. of this equation and 
it can easily be integrated. If we assume that 60 does not vary, i.e. 60 = — l, then 
the integral over dt is (approximately) proportional to (1 — q), so that for g < 1 the 
asymmetry decreases, and for g > 1 it starts to rise. It is easy to check that with 
D = the integrated rise is stronger than the decrease and thus the asymmetry rises 
by the factor exp [3.7(10^ sin 2^)^] for the mass Sm'^ = —1. For any sin2^ > 2 • 10~^ 
there would be an enormous rise of the asymmetry. 

This does not happen, however, because for a large 9 the variation of 60 should 
be taken into account. It changes as: 

Abo = &o(oo) - 6o(0) ^ 0.04 f^) \6m'\'/' , (54) 



Even this relatively weak variation happens to be vitally important for the evolution of 
the asymmetry. The point is that the integral in the r.h.s. of eq. ( |51| ) has a resonance 
at [/ = 0. The contribution of this resonance into the integral is quite small for 
a constant bo, because the factor U in the numerator cancels it out since U is an 
odd function near the resonance. However, since 69 experiences variation exactly at 
the resonance point, its variation breaks the symmetry from the positive and negative 
contribution of U near the resonance and the relative effect of the small variation of 69 
is enhanced by the factor 1/6 ~ 10^. This effect diminishes the positive contribution 
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into the r.h.s. of eq. (|5l|) and the rise of the asymmetry is strongly suppressed. In 
particular for sin 20 > 10^^ the integrated contribution of the r.h.s. becomes negative 
and the asymmetry decreases with respect to its initial value. 

The variation of 60 is also important because it gives an upper limit to a possible 
generation of lepton asymmetry. Because of leptonic charge conservation the asym- 
metry generated in the sector of active neutrinos must be equal to that in the sector 
of sterile ones. The latter is proportional to the difference of the diagonal matrix ele- 
ments of the density matrix, AZ ~ A(s — s). Since the variation of s and s can only 
be positive (initial value of both is —1), A(s — s) < A(s + s) and the last quantity is 
given by the variation of 60 • Naively taken from eq. (|5^ ) this variation is rather small 
for small 6. However, as we see in what follows, the variation of 60 is rising with rising 
Z, so the discussed limit is not broken. On the other hand, the back-reaction from 
the variation of bo terminates the rise of the asymmetry when it is still not too large; 
the maximum amplification, that happens to be near sin 20 ^ 10^^ for Sm^ = —1, 
could be about 10^. 

For 10^^ < sin 20 < 10^^ a solution of eq. (|5TD without back-reaction results in a 
huge rise of the asymmetry, however, the back-reaction efficiently kills this rise and 
the asymmetry may increase at most by 6 orders of magnitude. It confirms the early 
assertion (based on oversimplified arguments) of ref. that back-reaction does not 
permit the asymmetry to grow too much. However, the magnitude of the generated 
lepton asymmetry the we found is much larger than that advocated in ref. ^ but still 
much smaller than in refs. p[-||Tl|, except for a very large mass difference, 6m'^ ~ 10^ 
eV, where the asymmetry may be above 0.1. We have solved eq. (^TJ) numerically in 
the above mentioned range of sin 20. For sufficiently small values of q the equation 
was solved directly without any simplifications, while for larger q, when the product 
C = 3.3-10^^ Z{cos29Sm^)~^^^ became close to unity, the integral was estimated in the 
resonance approximation. There are two resonances corresponding to the condition 
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D"^ = U"^ that give opposite sign contributions to the integral, and eq. ( ^l]) becomes: 



c'= E ^&o(i/t 



^B, ,^,^,q''\t]-l)Uqt,) 



.=1,2 2 C\/4g2/3 + c 



'2 



2 1 N2 



l-5r 



(t -1 



6H-^ 



(55) 



where t^ = ±C/2q^l^ + V 1 + (C/2g^/3)2_ 

For (5m^ = — 1 we present the evolution of r^ as a function of the decreasing 
temperature T for 3 different mixing angles in fig. §, where the lepton asymmetry is 
V = \vb/4: — {riy^ — ny^)/n^\. The solid line is for sin 20 = 1 ■ 10~^, and it is clearly 
seen that the asymmetry is frozen at a very low value. For bigger mixing angles, 
sin26' = 2 ■ 10^^ (dashed) or sin26' = 3 ■ 10~^ (dotted) the increase may be much 
bigger. The dotted line clearly shows a power law behavior, t] oc T^^, following the 
exponential increase. If one neglects the back-reaction, i?i = 0, then the power law 
becomes r] oc T^^^^^. 

In fig. 1^ we plot the final value of r^ as a function of sin26' for 5m? = — 1 . One clearly 
sees the sharp exponential cut-off around sin20 = 10^^. The final lepton asymmetry 
in the region with a large increase, 10~^ < sin20 < 10^^, does not depend on the 
initial value, r/^, whereas the final value of t] is almost linear in rjin for sin26' < 10~^. 

For smaller masses the region of increase is shifted slightly to higher mixing an- 
gles as is seen in fig. ^, where we plot the final lepton asymmetry, r/, as functions of 
sin26' for 5 different masses, —5m? = 10^^, 1, 10^, 10^, 10^^. This shift is caused by 

,2hl/6 



B ~ (I (5m I) (see eq. (0)). On the other hand, for bigger masses the exponen- 



tial cut-off moves to smaller mixing angles. This is because A6o goes like J\5m,'^\. 
Clearly the effect with the biggest masses has limited applicability, since very heavy 
particles would become non-relativistic early, however, it is comforting to note that 
the maximal asymmetry generated does not continue rising for very heavy neutrinos. 
The region of instability in the {5m?, sin 26')-plane is presented in fig. ^. 

Our eq. ( ^Tf ) is very similar to the equations describing the evolution of the asym- 
metry derived in ref. |jT9|. However, we took terms related to the variation of bo 



into account and these terms are responsible for the stabilization of the rise of the 
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asymmetry when the latter is still small. Our derivation of the evolution equation is 
somewhat different and to our mind it is more rigorous. 

The result are only valid in the case of small mixing angles, sin 20 < 6 = 1/135. 
In the other limiting case the evolutionary equation is quite different as well as the 
behavior of the asymmetry. Our preliminary results show that in the case of a large 
mixing angle the asymmetry does not rise. This agrees with the result found in this 
work that for sin 26 > 10^^ the asymmetry diminishes. Even though the asymmetry 
remains small, the impact of neutrino oscillations on the primordial nucleosynthesis 
could be non-negligible due to the effects of non-equilibrium neutrino distribution 
function. We will calculate the abundances of light elements in a subsequent paper. 

Thus, according to our calculations, lepton asymmetry in the sector of active 
neutrinos may indeed be strongly enhanced. However the enhancement that we found 
is considerably weaker than that found in the earlier papers [^-[|TU|. Even for a very 



large mass difference —5m? = 10^ the resulting asymmetry is below 10^^ and only 
for —5ir? = 10^ it reaches the values that may be important for nucleosynthesis (see 
fig. ^. We believe that our calculations are very accurate. The only approximation is 
an expansion in inverse powers of Q (eq. (p9D) and in powers of a small sin 20 < 0.01. 
Otherwise our calculations are exact. In some cases we used analytic solutions to 
appropriate differential equations but in all the cases numerical solutions were quite 
simple and they agree very well with the analytic results. In the limit of a small 5m?, 
e.g. —5m? = 10^^ the parameter Q ~ 18 still remains large. Our results for such 
small 5m? are in a reasonable agreement with direct numerical calculations made in 
refs. |TB| and with our own ones (unpublished). 



We see that for some values of the mixing angle the asymmetry Z first very quickly 
goes to zero reaching extremely small values about or below lO^^'^'' and later started 
to rise up to 10~^ or even somewhat larger. If one solves eq. (|5l|) the sign of the 
final asymmetry remains fixed and is completely determined by the initial conditions. 
So in this approximation we do not observe any chaoticity in agreement with earlier 
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papers [^ |T4|- On the other hand if one apphes a direct numerical approach, then it is 
practically evident that chaoticity must be observed because no direct computational 
procedure is able to maintain an accuracy at the level 10^^°°. Thus one would expect 
that in the region when the asymmetry is quite small its numerically calculated sign 
is arbitrary and chaotic; it is just numerical errors. When the asymmetry starts to 
rise its final sign is the same as the initial one at the moment when the asymmetry 
becomes larger than the numerical accuracy. It could possibly explain the chaoticity 
observed in the papers [jTI], |T^, [T^ . 



However one should remember that eq. (|5T|) is valid only if Z does not decrease 
too fast with increasing q. So for a small Z one cannot say on the basis of eq. ( pT] ) 
that the asymmetry is as small as 10~^°°. There are some more terms in eqs. (p5|-p8|), 
as e.g. F(^A — S')/2, that should not be neglected if Z vanishes. Our preliminary 
results show that the solution of the kinetic equations in the limit of small Z does 
not show any chaoticity, though the sign of Z may be different from the initial one. 
We do not observe the sign change but this statement demands some further checks. 

There is a physically interesting possibility of chaoticity, namely if the asymmetry, 
as calculated through kinetic equations, is extremely small, the statistical fluctuations 
would be essential. The relative magnitude of a statistical fluctuation in a volume 
with A^ particles is about 1/yN . So if this value is larger than the asymmetry Z the 
fluctuations would dominate and the sign of the asymmetry would be determined by 
statistical fluctuations. However, to be essential the size of the region with such a fluc- 
tuation should be larger than the neutrino diffusion length during the characteristic 
time of oscillations. The complexity of the calculations in such a case would increase 
very much because now one has to take into account the effects of the fluctuating 
medium on the oscillations. 

A possible explanation of the difference between our results and the results of other 
groups (they also disagree between themselves, in particular in a possible chaotic be- 
havior of the asymmetry) is that in most cases an assumption of kinetic equilibrium 
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of neutrinos was made. This assumption enormously simplifies the numerical calcu- 
lations but may be strongly violated. Its violation may be crucial for the strength 
of generation of lepton asymmetry. We checked in a simplified example that in the 
opposite limit when the spectrum of neutrinos never recovered its equilibrium distri- 
bution and the resonance is complete (i.e. for relatively large sin 26'), the asymmetry 
experiences only a very mild enhancement fl% . 



However in some papers (see e.g. ||T0|) it is stated that a complete set of kinetic 
equations was numerically solved without any approximations. It is always difficult 
to find the source of disagreement, especially in numerical works. As we understood 



from the paper |T0[, the system of 8A^ kinetic equations (where N is the number of 



points in the momentum grid) for u and z/ (equivalent to our eqs. ([l^)-(lT7|)) was 
solved numerically but an additional equation for the evolution of lepton asymmetry 
was introduced {Z in our notations and L in notations of the quoted paper). The 
latter was obtained from the expression for the asymmetry by differentiating the cor- 
responding integrand containing elements of the density matrix. This equation was 
also solved numerically step by step and the resulting asymmetry was substituted 
into the equations describing the evolution of the density matrix elements. Possibly 
this technical trick helped to diminish the computational instability of the original 
equations. To calculate the integral the authors estimated it close to the resonance 
and integrated over the range of 3.5 resonance widths. We repeated a similar pro- 
cedure for the derivative of the asymmetry Z' in the resonance approximation and 
found both analytically and numerically that the result strongly depends upon the 
integration limits. The larger are the limits the weaker is the rise of asymmetry. The 
integral becomes saturated when the limits are larger than 10 or even 20 widths. Such 
a slow rate of saturation is related to the fact that the derivative of the resonance is a 
function that changes sign near the resonance. Another source of disagreement may 



be that the integration over momentum performed in ref. [^ was symmetric around 



the resonance, however in reality the range of integration in negative and positive 



directions are not the same (because the momentum runs from to oo). The term 
of the lowest order in the resonance width in the integrand is an odd function near 
the resonance so the contribution of asymmetry in the integration range is enhanced 
by 1/5. This effect is not strong in the case when the lepton asymmetry remains 
small and its back-reaction is not essential but when it starts to rise, the asymmetry 
in integration limits should be taken into account. These two effects may possibly 



explain the difference between our results and the calculations of ref. [|1^. However in 
our discussion with R. Foot and R. Volkas they defended stability of their calculations 
with respect to the choice of the region of integration. So at the moment the question 
about the precise origin of our disagreement remains open. 

In conclusion, we have analytically transformed the complete set of momentum 
dependent equations governing the evolution of the neutrino distribution functions 
to a form which allows a simple numerical solution. The only approximation is an 
expansion in the small parameter sin2^. These equations can even be solved analyt- 
ically in the limit of large Q, allowing us to derive a simple first order differential 
equation for the evolution of the lepton asymmetry. This differential equation takes 
into account the strong back-reaction effects on the generation of the lepton asym- 
metry due to the presence of an extra term (proportional to Bi) which is absent in 
approximate equations derived in some other papers. Due to this back-reaction we 
find that the asymmetry rise terminates at a much smaller magnitude. 
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Figure Captions: 

Fig. 1 bj as functions of q for momentum y = 1, sin2^ = 10^^ and Sm'^ = —1. 

The long-dashed hne is 1 + bo, the full line is 6i, the dotted and the dashed cures are 

absolute values of 62 and 63 respectively. 

Fig. 2 The evolution of r^ as a function of the decreasing temperature T in MeV. 

The full line is for sin2^ = 1 ■ 10~^, the dashed line is for sin2^ = 2 ■ 10~^, and the 

dotted line is for sin2^ = 3 ■ 10~^. All with Sm^ = —1. 

Fig. 3 The final value of r^ as a function of sin2^ for drri^ = — 1. For mixing angles 

bigger than ^ 10~^ the final value of r] is exponentially suppressed (see text). We 

have used (5m^ = —1, ?7j„ ~ 10"^*^ and q runs from 10^^ to 10^. 

Fig. 4 The final value of 77 as a function of sin2^ for 5 different masses: —5m? = 

10^6 (solid), 1 (dashed), 10^ (dotted), 10^ (dash-dot), and 10^^ (long-dashed). 

Fig. 5 Stability regions in sin 29 — 6ir? space. Here instability means that the 

final f] is more than an order of magnitude bigger than the initial 77. 
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